function [z, sc, Cs]=scoord_xy ( h, theta_s, theta_b, Tcline, N, z_type, index);
% SCOORD_XY:  computes depth of RHO- or W-points for a set of indices
%
% USAGE:  [z,sc,Cs]=scoord_xy(h,theta_s,theta_b,Tcline,N,z_type,index);
%
% PARAMETERS:
% Input:
%    h         Bottom depth (m) of RHO-points (matrix).                 
%    theta_s   S-coordinate surface control parameter (scalar):         
%                [0 < theta_s < 20].                                    
%    theta_b   S-coordinate bottom control parameter (scalar):          
%                [0 < theta_b < 1].                                    
%    Tcline    Width (m) of surface or bottom boundary layer in which  
%              higher vertical resolution is required during streching 
%              (scalar).                                               
%    N         Number of vertical levels (scalar).                     
%    z_type     Depth grid type logical switch:                         
%                z_type = 's_rho'   ->  depths of RHO-points.                
%                z_type = 's_w'   ->  depths of W-points.                  
%    index     Index into h matrix.
% Output:                                                              
%    z       Depths (m) of RHO- or W-points (matrix).                  
%    sc      S-coordinate independent variable, [-1 < sc < 0] at       
%            vertical RHO-points (vector).                            
%    Cs      Set of S-curves used to stretch the vertical coordinate  
%            lines that follow the topography at vertical RHO-points  
%            (vector).                                                
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% $Name: rslice-0_9_2 $
% $Id: scoord_xy.m,v 1.2 2005/05/26 21:46:09 jevans Exp $
% AUTHOR:  johnevans@acm.org
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 7
	help (mfilename);
	sprintf ( msg, '%s:  must have at least 7 input parameters.\n', mfilename );
	error ( msg );
end

%----------------------------------------------------------------------------
%  Set several parameters.
%----------------------------------------------------------------------------

c1=1.0;
c2=2.0;
p5=0.5;
Np=N+1;
ds=1.0/N;
hmin=min(min(h));
hmax=max(max(h));
hc=min(hmin,Tcline);
[Mp Lp]=size(h);

%----------------------------------------------------------------------------
% Define S-Curves at vertical RHO- or W-points (-1 < sc < 0).
%----------------------------------------------------------------------------

switch ( z_type )
case 's_rho'
	Nlev=N;
	lev=1:N;
	sc=-c1+(lev-p5).*ds;
case 's_w'
	Nlev=Np;
	lev=0:N;
	sc=-c1+lev.*ds;
otherwise
	msg = sprintf ( '%s:  unhandled z_type %s.\n', mfilename, z_type );
	error ( msg );
end


Ptheta=sinh(theta_s.*sc)./sinh(theta_s);
Rtheta=tanh(theta_s.*(sc+p5))./(c2*tanh(p5*theta_s))-p5;
Cs=(c1-theta_b).*Ptheta+theta_b.*Rtheta;
Cd_r=(c1-theta_b).*cosh(theta_s.*sc)./sinh(theta_s)+ ...
     theta_b./tanh(p5*theta_s)./(c2.*(cosh(theta_s.*(sc+p5))).^2);
Cd_r=Cd_r.*theta_s;

%============================================================================
% Compute depths at requested grid section.  Assume zero free-surface.
%============================================================================

zeta=zeros(size(h));

z=zeros(Nlev,length(index));
for k=1:Nlev,
	z(k,:)=zeta(index)'*(c1+sc(k))+hc.*sc(k) + (h(index)'-hc)*Cs(k);
end

z = flipud(z);


return


